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ABSTRACT We briefly review and highlight the conse- 
quences of rigorous and exact results obtained in [ ], char- 
acterizing the statistics of spike trains in a network of 
leaky Integrate-and-Fire neurons, where time is discrete 
and where neurons are subject to noise, without restriction 
on the synaptic weights connectivity. The main result is 
that spike trains statistics are characterized by a Gibbs dis- 
tribution, whose potential is explicitly computable. This es- 
tablishes, on one hand, a rigorous ground for the current in- 
vestigations attempting to characterize real spike trains data 
with Gibbs distributions, such as the Ising-like distribution 
[2], using the maximal entropy principle. However, it tran- 
spires from the present analysis that the Ising model might 
be a rather weak approximation. Indeed, the Gibbs poten- 
tial (the formal "Hamiltonian") is the log of the so-called 
"conditional intensity" (the probability that a neuron fires 
given the past of the whole network [3, 4, 5, 6, 7, 8, 9, 10]). 
But, in the present example, this probability has an infinite 
memory, and the corresponding process is non-Markovian 
(resp. the Gibbs potential has infinite range). Moreover, 
causality implies that the conditional intensity does not de- 
pend on the state of the neurons at the same time, ruling out 
the Ising model as a candidate for an exact characterization 
of spike trains statistics. However, Markovian approxima- 
tions can be proposed whose degree of approximation can 
be rigorously controlled. In this setting, Ising model ap- 
pears as the "next step" after the Bernoulli model (inde- 
pendent neurons) since it introduces spatial pairwise cor- 
relations, but not time correlations. The range of validity 
of this approximation is discussed together with possible 
approaches allowing to introduce time correlations, with 
algorithmic extensions. 

KEY WORDS Spike trains statistics, Gibbs measures, In- 
tegrate and Fire Models, chains with complete connections, 
Markov approximations, Ising model. 



1 Introduction. 

At first glance, the neuronal activity is manifested by the 
emission of action potentials ("spikes") constituting spike 
trains. Those spike trains are usually not exactly repro- 
ducible when repeating the same experiment, even with a 
very good control ensuring that the experimental conditions 
have not changed. Therefore, researchers are seeking sta- 
tistical regularities in order to provide an accurate model 



for spike train statistics, i.e. a probability distribution fit- 
ting the experimental data and/or matching what neurosci- 
entists believe models should predict. However, obtaining 
statistical models from experimental data remains a diffi- 
cult task. 

It appears simpler to characterize spike trains statis- 
tics in neural networks models where one controls exactly 
the neural network parameters, the number of involved neu- 
rons, the number of samples, and the duration of the exper- 
iment (with a possible mathematical extrapolation to infi- 
nite time). Especially, producing analytical (and when pos- 
sible, rigorous) results on those statistics provide clues to- 
ward resolving experimental questions and new algorithms 
for data treatments. Here, we propose a complete and rigor- 
ous characterization of spike train statistics for the discrete- 
time leaky Integrate-and-Fire model with noise and time- 
independent stimuli. This framework affords extrapolation 
to more realistic neural networks models such as general- 
ized Integrate-and-Fire [11, 12]. Our results hold for finite- 
sized networks, and all type of synaptic graph structure are 
allowed. Also, we are not constrained by an ad hoc choice 
of the initial conditions distribution of membrane poten- 
tials; instead we propose a procedure where this distribu- 
tion is selected by dynamics and is uniquely determined. 
More precisely, we show that spike train statistics are char- 
acterized by a (unique) invariant probability distribution 
(equilibrium state) whatever the model-parameters values, 
which satisfies a variational principle (the maximal entropy 
principle) and is a Gibbs distribution whose potential is ex- 
plicitly computed. This has several deep consequences dis- 
cussed in this paper. 

2 Model definition. 

We consider a discrete-time (continuous state) Integrate 
and Fire model, introduced in [13, 14], whose dynamics 
is given by: 



Vt(t+1) = jVi (1 - ZlViiW+^WijZlVjW+Ii+aBBiit) 
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(1) 

where i = 1 ... N is the neuron index, Vi the mem- 
brane potential of neuron i (this is a continuous variable), 
7 G [0, 1 [ is the leak rate, Z(x) = 1 whenever x > 9 and 
Z(x) = otherwise, where 6 is the firing threshold, Wij is 
the synaptic weight from neuron j to neuron i, Ii is a con- 
stant external current, <tb > 0, and Bi (t) is a noise, namely 



the Bj(t)'s are Gaussian i.i.d. random variable with mean 
zero and variance 1. The term (1 — Z[V[\) corresponds to 
the reset of membrane potential whenever neuron i fires 
while the term Y^jLi WijZ{Vj] is the post synaptic poten- 
tial generated by the pre-synaptic spikes. 



3 Spike train statistics. 

To each membrane potential value Vi(t) we associate a 
variable u>i(t) = Z(Vi(t)). The "spiking pattern" of the 
neural network at time t is the vector uj(t) = (oJi{t)) i=1 : 
it tells us which neurons are firing at time t, {tOi(t) = 1) 
and which neurons are not firing at time t (u>i(t) = 0). We 
denote by cj* the sequence or spike block u(s) . . . u>{t). A 
bi-infinite sequence ui = {w(i)}tt^oo of spiking patterns 
is called a "raster plot". It tells us which neurons are fir- 
ing at each time t S TL. In practical experimental raster 
plots are obviously finite sequences of spiking pattern but 
the extension to Z, especially the possibility of consider- 
ing an arbitrary distant past (negative times) is a key tool 
in the present work. We are seeking invariant probability 
distributions on the set of raster plots, generated by the dy- 
namical system (1). In the next sections we give our main 
results whose proofs can be found in [1]. 

4 The transition probability. 

In this model it is possible to compute exactly and rigor- 
ously the probability, called conditional intensity in [3, 4, 
5, 6, 7, 8, 9, 10], of a spiking pattern at time t+ 1, uj(t+ 1), 
given the past of the network a/.^. To our knowledge this 
is the first time that such an exact and rigorous computation 
is achieved at the level of network. It is given by: 



N 



p^t+i)!^,) =n^M*+i)i 



(2) 



2 = 1 



with 



P(uji(t+ 1)1^) = 

u>i(t + l)ir{ ai{ut ] j+ (3 ) 
(l-^ + l))(l-,(^g^)), 



where 
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C^-oo) = £ WySii(<"-oo) + li ■ 

with 



'7 



^•j(^-oo) = J2 7 ' WO. 

is the deterministic contribution of the network to the mem- 
brane potential of neuron i at time t + 1, when the spike 



train up to time t is iJL x . The term Ci{uj t _ x ) integrates 
the pre-synaptic flux that arrived at neuron i in the past, as 
well as the external current Ij. As a consequence, this term, 
depends on the past history up to a time r^ui^) which is 
the last time before t when neuron i has fired, in the se- 
quence (the mathematical definition is given below). 



Likewise, 



Oi(W-ao) 
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is the variance of the noise-induced term applied to neu- 
ron i at time t + 1 and resulting from the integration of 
the noise from time t^cj^^) to time t. Finally, tt(x) = 



r 



e 2 dw. 



The mathematical definition of r^wl^) is: 



clef 



-oo, if w i (/c)=0, Vfc<t; 
max{-oo < k < t,uji(k) = 1} otherwise. 



Basically, this notion integrates the fact that the state of 
neuron i depends on spikes emitted in the past by the net- 
work, up to the last time when this neuron has fired (thus, 
the present analysis relies heavily on the structure of IF 
models where reset of the membrane potential to a fixed 
value erases the past evolution). Depending on the se- 
quence uLoo, this time can go arbitrary far in the past. 
Now, here, sequences such that r^wi^) < — c for some 
positive c, have a positive probability (decaying exponen- 
tially fast with c). The consequence is that the transition 
probability (2) depends on the history of the network on 
an unbounded past. As a consequence dynamics is non 
Markovian. 

Note also that (evidently ?) the structure of the neural 
network dynamics imposes causality. This appears explic- 
itly in (3) where the probability that neuron i fires at time 
t+1 depends on spikes emitted up to time t. In other words, 
this probability does not depend on the spikes emitted by 
the other neurons at the same time (t + 1). We shall come 
back to this remark when discussing the relevance of Ising 
model. 

Finally, remark that the factorization of the probabil- 
ity (2) expresses the conditional independence of spikes 
at given time. This is simply due to the fact that the only 
source of stochasticity, when past spikes are fixed (by 
the conditioning), is the noise, which is assumed, in this 
model, to be independent between neurons. But, (2) does 
not imply that spikes are independent. 

Let us emphasize what we have obtained. We have an 
explicit form for the probability that a neuron fires given 
its past. From this we can obtain the probability of any 
spike block. But, here there is technical difficulty, since the 
corresponding stochastic process is not Markovian. Such 
process corresponds to a more elaborated concept called a 
chain with compete connections. As a consequence the ex- 
istence and uniqueness of an invariant measure is not given 
by classical theorems on Markov chains but requires more 



elaborated tools, widely developed by the community of 
stochastic processes and ergodic theory (see [15] for a nice 
introduction and useful references). The price to pay is a 
rather abstract mathematical formalism but the reward is 
the possibility of considering spike statistics arising in such 
models, including memory effects and causality, with, we 
believe, possible extensions toward more realistic neural 
models. 

5 Gibbs distribution. 

In the present setting where Z, does not depend on t, P 
is stationary. Namely, fix a spiking sequence {a n } n<0 , 
then, Vi, P {u){t) = a \ u(t — ri) = a_ n , n > 1) = 
P (clq I aZoo) ■ Therefore, instead of considering a family 
of transition probabilities depending on t, it suffices to 
define the transition probability at one time t 6 7L, for 
example t = 0. The next results are based on techniques 
from ergodic theory and chains with complete connections 
[15]. 

We recall first a few definitions. A Gibbs distribution 
is a probability distribution fi^, on the set of (infinite) spike 
sequences (raster plots) X, such that one can find some 
constants P(ip), C\, c 2 with < ci < 1 < c 2 such that for 
all n > and for all u> € X: 



ci < 



W W ) 



exp[-(n+l)P(V)+ELo^( T M] " 



f'2- 



where T is the right shift 1 over the set of infinite sequences 
X i.e. (Tcj)(t) = uj(t - 1). P(ip), called the topologi- 
cal pressure, is formally analogous to the free energy and 
is a central quantity since it allows the computation of av- 
erages with respect to fj,^. In the previous equation, the 
term Efc=o V'C* may be viewed as a formal Hamilto- 
nian over a finite chain of symbols u!q. However, as exem- 
plified below (eq. (5)) this term depends also on the past 
(infinite chain ujZ^)- There is an analogy with statistical 
physics systems with boundary conditions where the Gibbs 
distribution not only depends on the configuration inside 
the lattice, but also on the boundary term [16]. 

Moreover, fj,^ is an equilibrium state if it satisfies 
("maximum entropy principle"): 

P(tp) d = h(m) + MV>(V>) = sup h(fx) + fj,(tp), (4) 

where Pr(X) the set of T-invariant finite measures on X, 

is the entropy of p,^, and /i(V>) = / ipd/i is the av- 
erage value of ip with respect to the measure fi. Though 
Gibbs states and equilibrium states are non equivalent no- 
tions in the general case, they are equivalent in the present 
setting. 

Now, our main result is: 



The dynamical system (1) has a unique invariant- 
measure, pL^, whatever the values of parameters 



Wij, i,j 
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1 . . . N, 7, 6. This is a Gibbs 



distribution and an equilibrium state, for the potential: 



#*)=log(P(w(0)|wli,))s^(a;° 
= Ef =1 h(0)log(vr(^M" 

(l-^(0))log(l-7T 



0-4(0;) 
8-Ci(iiJ) 



(5) 



where we note CJ 
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The use of the right shift instead of the left shift is standard in the 
context of chains with complete connections [15]. 



Note therefore that the Gibbs potential is simply the 
log of the conditional intensity. We believe that this last 
remark extends to more general models. 

Knowing that is Gibbs distribution with a known 
potential allows one to control the probability distribu- 
tion of finite sequences. Note that here the potential i/j 
is explicitly known. This has several other deep conse- 
quences. First, this formalism opens up the possibility of 
characterizing P(ip) and ^ by a spectral approach, be- 
ing respectively the (unique) largest eigenvalue and related 
left eigenfunction of the Ruelle-Perron-Frobenius operator 
AJH = E W ' : T(o/)=" acting on C(X, R), 

the set of continuous real functions on X. This last property 
has deep consequences at the implementation level [17]. 

Also, it is possible to propose Markovian approxima- 
tions of this distribution, whose degree of accuracy can be 
controlled, as we now discuss. 

6 Markov approximations. 

The main difficulty in handling the transition probabilities 
(2) and the related equilibrium state is that they depend on 
an history dating back to r^wL^), where ^(oA^) is un- 
bounded. On the other hand, the influence of the activ- 
ity of the network, say at time — I, on the membrane po- 
tential Vi at time 0, appearing in the term .t^^w^^) = 
Ei=7~fui° )l~ luj j(l)> decays exponentially fast as I — > 
— 00. Thus, one may argue that after a characteristic time 
depending on i \ g^\ the past network activity has little 
influence on the current membrane potential value. 

Assume that we want to approximate the statistics of 
spikes, given by the dynamics (1), by fixing a finite time 
horizon R such that the membrane potential at time de- 
pends on the past only up to some finite time — R. In 
this way, we truncate the histories and we approximate the 
transition probabilities P (w(0) |wZ^o), with unbounded 
memory, by transition probabilities P (w(0) | wZ/j), thus 
limiting memory to at most R time steps in the past. These 
approximated transition probabilities constitute therefore a 
Markov chain with a memory depth R. But how good is 
this approximation ? 

The truncation of histories leads to a truncation of the 
Gibbs potential, denoted The invariant measure of 

the Markov chain is a Gibbs distribution for the potential 



ip( R \ denoted fi^R). One can show that the Kullback- 
Leibler divergence between the exact measure and the 
Markov measure obeys: 

d (h^r) , < C-f R , (6) 

where C can be computed explicitly as a function of synap- 
tic weights and current. 

Therefore, the Kullback-Leibler divergence decays 
exponentially fast, with a decay rate 7, as expected from 
our prior qualitative analysis. It is thus sufficient, for prac- 
tical purposes, to approximate ip with a potential of range: 

R ~ it— — 1 • (7) 

I log Tl 

Consequently, the smaller the leak, the shorter is the range. 

Note that the previous result is classical in ergodic 
theory and expresses that finite range potentials are dense 
in the set of Gibbs potentials [16]. What we bring is the 
computation of the decay rate and an estimation of the con- 
stant C for the present model (see [ 1 ] for details). 

7 Raster plots statistics 

As discussed in the introduction, the neuroscience commu- 
nity is confronted to the delicate problem of characterizing 
statistical properties of raster plots from finite time spike 
trains and/or from finite number of experiments. This re- 
quires an a priori guess for the probability of raster plots, 
what we call a statistical model. These models can be ex- 
trapolated from methods such as Jaynes' [18]: "Maximis- 
ing the statistical entropy under constraints". Actually, this 
is none other that the variational formulation (4). Now, this 
method only provides an approximation of the sought prob- 
ability, relying heavily on the choice of constraints usually 
fixed from phenomenological arguments, and this approxi- 
mation can be rather bad [19]. Moreover, note that Jaynes' 
method is typically used for a finite number of constraints 
while (4) holds for a very general potential, actually cor- 
responding, in the present context, to infinitely many con- 
straints. 

On phenomenological grounds, spike n-uplets 
(ii, fi), . . . , (i n , t n ) (neuron i\ fires at time t\, and 
neuron 12 fires at time t^, ■ ■ ■ ) provide natural choices 
of constraints since they correspond to experimentally 
measurable events whose probability of occurrence has a 
phenomenological relevance. For example, the probability 
of occurrence of (ii,t\), /^(cjjl(ti)), is the instantaneous 
firing rate of neuron i\ at time t\. On a more formal ground 
a spike n-uplet corresponds to an order-n monomial; this 
is a function <j> which associate to a raster cj the product 
a;jj(ti) . . .u!i n (t n ), where 1 < i\ < i<x < • • • < i n < N 
and —00 < ti < ti < ■ ■ • < t n < 0, and where there is 
no repeated pair of indexes (i, t). Hence, </>(o;) is equal to 
1 if and only if neuron i\ fires at time t\, and neuron 12 
fires at time t^, . . . in the raster u> (and is otherwise). A 
polynomial is a linear combination of monomials. 



Fixing constraint on spike-uplets (monomials) aver- 
age leads to specific forms of Gibbs distribution. As an 
example, constraining firing rates only corresponds to a 
Bernoulli distribution where neurons spike independently 
and where the potential reads J^iLi ^i w ;(0)- Constrain- 
ing rates and probability that 2 neurons fire simultaneously 
leads to the so-called Ising distribution where the potential 
reads J2i=i ^ w i(0) + Z) l .j=i XijUi(0)cjj(0) [2]. Here, the 
related probability measure does not factorize any more but 
all information about spike train statistics is contained in 
the first and second order spike-uplets at time 0. 

More generally, imposing constraint on a (finite) 
set of monomials <f>\ , . . . , <f>L, leads to a parametric form 
of Gibbs potential (p guess (uj) = Ylh=i ^<M W )' where 
the A; are Lagrange multipliers. Note that the Az's are 
related by an additional constraint, the normalisation of 
the probability. Again, this potential relies heavily on the 
choice of constraints. 

On the opposite, in the present example, instead 
of an ad hoc guess, an exact polynomial expansion 
can be obtained from the explicit form of the potential 
(5). It is given by expanding (5) in series via the se- 
ries expansion of log(7r( 6 '~ c |'J^- ) )) (note that 8 ~J^^ > 
is bounded in absolute value provided the synaptic 
weights are finite and the noise intensity <tb is posi- 
tive. Thus, < a < ir( 9 ~~ C 7^ ) < b < 1, for some 
a,b). This series involves terms of the form C,(w) n 
leading to monomials of the form lj j 1 (ti) . . . u>j n (t n ) 
where ji , . . . , j n are neurons pre-synaptic to i (with 
W ijk ^ 0, . . . k = 1 . . . n) and t k < 0,k = l...n. 
As a consequence, the potential (5) has the form 
Yli=i Si^i ^i,i ( t ) i,i{ ul ) where the <fti/s are monomials 
of the form ^(0)^ (ti) . . . ojj n (t n ), t\,...t n < 
and the Au's contains products of the form 
^ ...^W ih ...W ijn . 

Note that terms of the form u>i 1 (0)wi 2 (0) . . . do not 
appear in this expansion. Especially, it does not contain the 
Ising term. 

As discussed in the previous section, truncating this 
series to a finite polynomial involving monomials with a 
time depth R (i.e. of the form ujj 1 {ti) . . .ujj n (t n ) with 
— 1 < tk < —R) amounts to considering Markovian ap- 
proximations of ip. How the corresponding Gibbs distri- 
bution approximates the exact one is controlled by eq. (6) 
where C depends on the synaptic weights and currents. 

This expansion contains an exponentially increasing 
number of terms as R, or N, the number of neurons, 
growths. However, it does not contain all possible mono- 
mials. Indeed, beyond, the remark above about the van- 
ishing of non causal term the form of the potential can be 
considerably reduced from elementary considerations such 
as stationarity. Moreover, some A^ ti i t( 's can be quite 
close to zero and eliminating them does not increase signif- 
icantly the Kullback-Leibler divergence. This can be used 
to perform systematic (algorithmic) reduction of the poten- 



tial in a more general context than the present model (see 
next section and [ 1 7] for more details). 

To summarize, one can obtain, from the previous 
analysis, a canonical form for a range-i? approximation of 
tfj, of the form: 

K 
2=0 

where all terms contribute significantly to the probability 
distribution i.e. removing them leads to a significant in- 
crease of the KL divergence between the exact Gibbs dis- 
tribution and its approximation. Thus, in the present case, 
the relevant constraints (and the value of the related A; 's) 
can be estimated from the analytic form of the potential. 

8 Parametric estimation of spike trains 
statistics. 

This analysis opens up the possibility of developing effi- 
cients algorithms to estimate at best the statistic of spike 
trains from experimental data, using several guess poten- 
tial and selecting the one which minimizes the KL diver- 
gence between the Gibbs measure [i , (rj and the empirical 

measure attached to some experimental raster oj^p) [17]. 
The idea is to start from a parametric form of potential (8), 
of range R, and to compute the empirical average of all 
monomials <pi from the experimental raster u)\- exp >. Then, 
one adjust the parameters A; by minimizing the KL diver- 
gence. This computation can be easily done using spectral 
properties of the Perron-Frobenius operator. This algorithm 
described in [17], will be presented in this conference, in 
another communication, and is freely available as a C++ 
source at http : / / enas . gf orge . inria . f r. 

9 Discussion 

In this paper we have addressed the question of character- 
izing the spike trains statistics of a network of LIF neurons 
with noise, in the stationary case, with two aims. Firstly, 
to obtain analytic and rigorous results allowing the charac- 
terization of the process of spike generations. Secondly, to 
make a connection from this mathematical analysis toward 
the empirical methods used in neuroscience community for 
the analysis of spike trains. Especially, we have shown that 
the "Ising model" provides, in this example, a bad approx- 
imation of the exact Gibbs potential. This is due, on one 
hand, to the fact that Ising potential does not take into ac- 
count causality, and on the other hand, to the fact that the 
exact potential includes quite a bit more "constraints" than 
the mere average value of pairwise spike coincidence. Al- 
though the first objection can somewhat be relaxed when 
considering data binning (that we do not consider here), 
the second one seems unavoidable. 

Also, we have shown that the Jaynes method, based 
on an a priori choice of a "guess" potential, with finite 
range, amounts to approximate the exact probability 



distribution by the Gibbs distribution of a Markov chain. 
The degree of approximation can be controlled by the 
Kullback-Leibler divergence. 

One may wonder whether our conclusions, obtained 
for a rather trivial model from the biological point of view 
has any relevance for the analysis of real neurons statistics. 
We would like to point that the main ingredients making 
this model-statistics so complex are causality induced by 
dynamics, and integration over past events via the leak 
term. We don't see any reason why the situation should 
be simpler for more realistic models. Especially, what 
makes the analysis tractable is the reset of the membrane 
potential to a constant value after firing, inherent to IF 
models, and rather unrealistic in real neurons. Thus, the 
memory effects could be even worse in realistic neurons, 
with a difficulty to extract, from a thorough mathematical 
analysis, the relevant times scales for a memory cut-off, 
as the log of the leak is in the present model. Neverteless, 
this work is a beginning with a clear need to be extended 
toward more realistic models. 

A natural extension concerns the so-called General- 
ized Integrate and Fire models [11] , which are closer to 
biology [21]. The occurrence of a post-synaptic potential 

(n) 

on synapse j, at time tj. , results in a change of membrane 
potential. In conductance based models this change is in- 
tegrated in the adaptation of conductances. It has been 
shown in [12], under natural assumptions on spike-time 
precision, that the continuous-time evolution of these equa- 
tions reduces to the discrete time dynamics. In this case the 
computation of the potential corresponding to (5) is clearly 
more complex but still manageable. This case is under cur- 
rent investigations. 

One weakness of the present work is that it only con- 
siders stationary dynamics, where e.g. the external current 
Ii is independent of time. Besides, we consider here an 
asymptotic dynamics. However, real neural systems are 
submitted to non static stimuli, and transients play a crucial 
role. To extend the present analysis to these case one needs 
the proper mathematical framework. The non stationar- 
ity requires to handle time dependent Gibbs measures. In 
the realm of ergodic theory applied to non equilibrium sta- 
tistical physics, Ruelle has introduced the notion of time- 
dependent SRB measure [ ]. A similar approach could be 
used here, at least formally. 

In neural networks, synaptic weights are not fixed, 
as in (1), but they evolve with the activity of the pre- and 
post-synaptic neuron (synaptic plasticity). This means that 
synaptic weights evolve according to spike train statistics, 
while spike train statistics is constrained by synaptic 
weights. This interwoven evolution has been considered 
in [23] under the assumption that spike-train statistics is 
characterized by a Gibbs distribution. It is shown that 
synaptic mechanism occurring on a time scale which is 
slow compared to neural dynamics are associated with a 
variational principle. There is a function, closely related 



to the topological pressure, which decreases when the 
synaptic adaptation process takes place. Moreover, the 
synaptic adaptation has the effect of reinforcing specific 
terms in the potential, directly related to the form of the 
synaptic plasticity mechanism. The interest of this result 
is that it provides an a priori guess of the relevant terms in 
the potential expansion. A contrario, it allows to constrain 
the spike train statistics of a LIF model, using synaptic 
plasticity with an appropriate rule which can be determined 
from the form of the expected potential. 

Acknowledgment:. We are indebted to T. Vieville and H. 
Rostro for valuable discussions and comments. 
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